Journal  of  Power  Sources  183  (2008)  98-108 


ELSEVIER 


Contents  lists  available  at  ScienceDirect 

Journal  of  Power  Sources 

journal  homepage:  www.elsevier.com/locate/jpowsour 


pri 

SbbaautS 


Control-orientated  thermal  model  for  proton-exchange  membrane 
fuel  cell  systems 

G.  Vasu,  A.K.  Tangirala* 

Indian  Institute  of  Technology-Madras,  Chennai  600036,  India 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  11  December  2007 
Accepted  25  March  2008 
Available  online  18  April  2008 


Keywords: 

Proton-exchange  membrane  fuel  cell 
Water  pump 
Thermal  model 
Interaction  analysis 
Stack  temperature 


A  lumped  parameter  dynamic  model  is  developed  for  predicting  the  stack  temperature,  temperatures  of 
the  exit  reactant  gases  and  coolant  water  outlet  in  a  proton-exchange  membrane  fuel  cell  (PEMFC)  system. 
A  dynamic  model  for  a  water  pump  is  also  developed  and  can  be  used  along  with  the  thermal  model  to 
control  the  stack  temperature.  The  thermal  and  water  pump  models  are  integrated  with  the  air  flow 
compressor  and  PEMFC  stack  current-voltage  models  developed  by  Pukrushpan  et  al.  to  study  the  fuel 
cell  system  under  open  and  closed-loop  conditions.  The  results  obtained  for  the  aforementioned  variables 
from  open-loop  simulation  studies  are  found  to  be  similar  to  the  experimental  values  reported  in  the 
literature.  Closed-loop  simulations  using  the  model  are  carried  out  to  study  the  effect  of  stack  temperature 
on  settling  times  of  other  variables  such  as  stack  voltage,  air  flow  rate,  oxygen  excess  ratio  and  net  power 
of  the  stack.  Further,  interaction  studies  are  performed  for  selecting  appropriate  input-output  pairs  for 
control  purpose.  Finally,  the  developed  thermal  model  can  assist  the  designer  in  choosing  the  required 
number  of  cooling  plates  to  minimize  the  difference  between  the  cooling  water  outlet  temperature  and 
stack  temperature. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cells  are  electrochemical  reactors  which  directly  con¬ 
vert  chemical  energy  to  electrical  energy.  Fuel  cells  are  typically 
classified  based  on  the  type  of  electrolyte  used  in  the  cell.  The 
proton-exchange  membrane  fuel  cell  (PEMFC)  is  widely  studied 
and  attains  its  name  from  the  use  of  a  polymer  membrane  as  the 
electrolyte.  Most  of  the  current  research  and  development  activ¬ 
ity  focuses  on  the  PEMFC  due  to  its  versatile  features  such  as  high 
power  density,  relatively  fast  start-up,  and  short  response  times  to 
changes  in  the  power  demand. 

In  automobile  applications,  an  important  requirement  is  that  the 
stack  should  meet  the  load  demands  of  a  varying  profile  with  short 
transient  times.  In  this  regard,  knowledge  of  the  stack  in  terms  of 
steady  state  and  transient  behaviour  is  of  critical  importance.  The 
dynamic  behaviour  of  a  PEMFC  system  is  strongly  dependent  on 
the  reactant  flows  and  the  water  and  thermal  management.  One  of 
the  several  challenges  that  arise  in  the  control  of  a  PEMFC  system 
is  the  level  of  interaction  among  these  various  factors  such  that  an 
interaction  study  is  necessary  [1  ]  to  understand  the  relative  impor¬ 
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tance  of  these  parameters.  In  fact,  water  and  thermal  management 
have  become  the  key  technical  challenges  for  fuel  cell  technology 
to  be  feasible  for  transportation  applications.  Proper  thermal  and 
water  management  is  in  fact,  also  essential  to  achieve  optimum 
performance  from  PEMFC  stacks  [2,3]. 

Mathematical  models  play  an  important  role  in  supporting  the 
design  and  enhancing  the  understanding  of  the  effect  of  parame¬ 
ters  on  the  performance  of  the  stack  and  fuel  cell  auxiliary  systems. 
Studies  based  on  these  models  are  useful  for  the  optimum  design 
and  control  of  a  real-time  stack.  Two  modeling  approaches  can  be 
found  in  the  literature.  The  first  gives  rise  to  what  are  known  as 
mechanistic  models  [4-6]  which  use  an  in-depth  knowledge  of  the 
electrochemistry,  heat  transfer  and  mass  transfer  that  are  involved 
in  the  fuel  cells.  Such  models  explain  the  fundamental  processes 
occurring  in  fuel  cell  systems,  and  are  developed  as  ID,  2D  and 
3D  models  depending  on  the  assumptions  involved  therein.  These 
dimensional  models  for  thermal  and  water  management,  which  are 
summarized  in  [7,8],  require  iterative  methods  to  solve  the  under¬ 
lying  differential  and  partial  differential  equations,  thereby  making 
them  computationally  intensive.  In  essence,  the  mechanistic  mod¬ 
els  are  suited  for  design  and  optimization  of  the  individual  cell 
components,  rather  than  for  control  and  monitoring  of  the  stack. 
The  second  approach  includes  models  that  are  based  on  empirical 
or  semi-empirical  equations  which  are  applied  to  predict  the  effect 
of  different  input  parameters  on  the  voltage-current  characteris¬ 
tics  of  the  fuel  cell  [9-11].  When  compared  with  the  mechanistic 
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Nomenclature 

A  area(m2) 

Cp  specific  heat  capacity  (J  mol-1  K-1 ) 

g  gravity  constant  (ms-2) 

h  heat  transfer  coefficient  (W  m-2  K-1 ) 

AH  heat  of  reaction  (J  mol-1 ) 

I  current  (A) 

J  moment  of  inertia  (kg  m2 ) 

k  thermal  conductivity  of  air  (W  m-1  K-1 ) 

L  length  of  stack  (m) 

m  mass  (kg) 

M  molar  mass  (mol  g-1) 

n  number  of  cells  in  stack 

N  molar  flow  rate  (mol  s-1) 

Npr  Prandtl  number 

PEM  proton-electrolyte  membrane 

Q.  water  flow  rate  (L s-1 ) 

Q  energy  (W) 

Ra  Rayleigh  number 

T  temperature  (K) 

V  stack  voltage  (V) 

W  mass  flow  rate  (kg  s-1 ) 

Greek  letters 

v  kinematic  viscosity  (m2  s-1 ) 

p  density  (kg  m-3) 

o)r  motor  speed  (rps) 

Subscripts 
a  anode 

c  cathode 

cell  proton-exchange  membrane  cell 

com  combined 

con  consumed 

elec  electrical  energy 

H2  hydrogen 

H20  Water 

in  inlet 

1  liquid 

latent  latent  heat 

LMTD  Logarithmic  mean  temperature  difference 

max  maximum 

mot  motor-pump 

N2  nitrogen 

out  outlet 

02  oxygen 

sens  sensible  heat 

stack  proton-exchange  membrane  fuel  cell  stack 

theo  theoretical  chemical  energy 

trans  net  transfer  of  water 

v  vapour 

W  water 

0  standard  condition 


models  the  physical  and  electrochemical  phenomena  are  modelled 
at  coarser  levels.  The  semi-empirical  models  for  thermal  manage¬ 
ment  of  a  stack  cannot  be  used  directly  for  control  studies  due 
to  the  implicit  forms  of  the  expressions  that  are  used  to  calculate 
the  stack  temperature.  From  a  control  perspective,  therefore,  it  is 
important  to  develop  a  model  that  describes  the  thermal  behaviour 
of  a  PEMFC  stack  with  good  approximation  of  the  dynamics  of  stack 
temperature  under  various  load  conditions. 


The  total  thermal  energy  evolved  during  operation  of  a  fuel  cell 
is  the  difference  between  the  chemical  energy  of  H2-02 (input)  and 
the  electrical  energy  of  the  stack  (output).  This  thermal  energy  is 
distributed  as  sensible  heat  of  the  coolant  and  reactants,  latent  heat 
during  phase  change  of  water  and  heat  loss  to  the  surroundings  by 
convection.  Most  of  the  lumped  parameter  models  that  have  been 
developed  for  predicting  the  stack  temperature  have  considered  the 
sensible  heat  of  cooling  water  as  a  sole  transport  factor  of  the  ther¬ 
mal  energy  [12,13].  One  such  simple  thermal  model  of  the  PEMFC 
stack  was  embedded  into  the  ADVISOR  vehicle  simulation  pack¬ 
age  by  group  of  researchers  at  NREL  [14].  The  required  inputs  and 
parameters  for  simulating  that  model  were  taken  from  the  ADVISOR 
package  of  automotive  fuel  cell  stack  driving  cycles  [15]. 

A  paradigm  shift  is  observed  in  the  work  by  Yu  et  al.  [16]  where 
the  latent  heat  of  water  during  the  phase  change  of  water  in  the  fuel 
cell,  the  heat  loss  to  the  surroundings,  the  sensible  heat  of  coolant 
water  and  reactants  are  taken  into  account  to  predict  the  temper¬ 
atures  of  the  stack,  the  exit  reactant  gases  and  the  coolant  outlet 
[  16].  In  this  model,  the  heat  loss  to  the  surroundings  by  convection  is 
not  considered  explicitly  but  rather  is  modelled  in  an  implicit  way. 
This  implicit  form  again  limits  the  utility  of  the  model  for  control 
application. 

In  the  present  work,  control-orientated  system-level  dynamic 
models  are  developed  for  (i)  stack  temperature  dynamics  by  explic¬ 
itly  taking  into  account  the  heat  loss  to  the  surroundings  in  addition 
to  the  latent  heat  of  vapourization,  the  sensible  heat  of  coolant 
water  and  reactants,  and  (ii)  the  centrifugal  water  pump  that  is 
used  to  control  the  stack  temperature.  These  models  are  inte¬ 
grated  with  first-principle  based  models  of  an  air  flow  compressor 
and  the  semi-empirical  voltage-current  model  of  Pukrushpan  et  al. 
[17]. 

The  results  of  the  present  work  offer  design  ideas  to  (i)  achieve 
the  exit  coolant  water  temperature  and  near  stack  temperature 
and  (ii)  decrease  the  settling  time  of  the  stack  temperature.  To 
the  best  of  the  authors’  knowledge  there  has  been  no  work 
reported  in  the  literature  that  addresses  the  development  of 
a  control-orientated  thermal  model.  Besides,  the  water  pump 
sub-system  model  provides  a  realistic  setting  for  studying  the 
temperature  control  problem.  The  goals  of  the  present  work 
are  to 

•  develop  a  control-orientated  lumped  parametric  system  level 
model  for  the  stack  temperature; 

•  develop  a  lumped  parametric  model  for  the  water  pump  sub¬ 
system; 

•  design  the  number  of  cooling  plates  required  to  achieve  the  cool¬ 
ing  water  exit  temperature  close  to  the  stack  temperature; 

•  study  the  interaction  between  inputs  (compressor  motor  voltage, 
motor-pump  current)  and  outputs  (air  flow  rate,  stack  tem¬ 
perature)  using  developed  thermal  and  water  pump  models  in 
conjunction  with  stack  and  air  flow  sub-system  models; 

•  study  the  effect  of  stack  temperature  on  the  settling  times  of  air 
flow  rate,  stack  voltage  and  net  power. 

This  paper  is  organized  as  follows.  First,  the  development  of 
a  dynamic  model  of  the  centrifugal  water  pump  is  explained  in 
Section  2.  The  proposed  control-orientated  thermal  management 
model  is  then  elucidated  in  Section  3.  A  brief  review  of  the  air  flow 
compressor  and  the  PEMFC  stack  models,  which  are  used  in  build¬ 
ing  the  integrated  model,  is  provided  in  Appendix  A.  Subsequently, 
in  Section  4,  the  results  from  simulation  studies  are  presented  with 
insights  into  the  number  of  cooling  plates,  the  knowledge  of  which 
is  required  during  the  design  of  the  stack.  Concluding  remarks  and 
directions  for  future  work  are  given  in  Section  5. 
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2.  Development  of  water  pump  model 

A  dynamic  lumped  parameter  model  is  developed  for  the  water 
pump  with  the  objective  of  controlling  the  stack  temperature  by 
manipulating  the  coolant  water  flow  rate.  The  coolant  water  flow 
rate  can  be  manipulated  in  two  different  ways,  namely:  (i)  by  using 
a  control  valve  in  combination  with  a  constant  speed  pump;  (ii)  by 
using  a  variable  speed  pump.  The  second  arrangement  is  attractive 
since  it  does  not  involve  the  use  of  a  valve  and  hence  reduces  the 
complexity  and  the  weight  of  the  system.  In  view  of  this,  dynamic 
modelling  of  the  variable  speed  centrifugal  pump  is  taken  up.  The 
water  flow  rate  of  the  pump  is  modelled  as  a  function  of  motor- 
armature  current  and  motor  speed.  The  water  pump  model  along 
with  the  thermal  model  facilitates  a  rigorous  and  realistic  study  of 
the  stack  temperature  control  problem.  It  is  worth  noting  that  the 
dynamics  of  the  water  pump  have  been  neglected  in  earlier  investi¬ 
gations  pertaining  to  control  of  the  stack  temperature  [12,13,16,18]. 

The  dynamic  pump  model  is  developed  based  on  the  fundamen¬ 
tal  relationships  between  motor-armature  current  (input),  motor 
speed  and  water  flow  rate  (Q)  (output)  of  the  centrifugal  pump 
reported  by  Kallesge  [19].  The  equations  governing  the  relation- 


ships  among  the  variables  are 

do)r 

Jcom  ^  —  Te  —  dCOt  -  ip 

(i) 

JwW=Hp-Pl 

(2) 

Hp  =  — fl|Q.3  +  q2Qo)y  +  ci3co2 

(3) 

Tp  =  —  fq  Q.3  +  b2Qoi>r  +  b3co2 

(4) 

P\  =  Pout  —  P[n  +  psiz  out  —  Z[n)  —  {Kv  +  Kp)Q2 

(5) 

where  Pi ,  P2,  Zi ,  z2  are  the  pressure  and  height  values  at  the  surface 
of  the  reservoir  and  the  stack  inlet,  respectively. 

The  combined  torque  (re)  equation  for  a  dc  motor  with  varying 
armature  current  is  given  by  [20] 

Te  =  C/mot  —  Tf  (6) 

The  equation  for  load  pressure  provided  above  in  Eq.  (5)  is  mod¬ 
ified  for  a  centrifugal  pump  of  constant  head  based  on  the  following 
assumptions: 

•  No  control  valve  is  used  to  control  the  water  flow  rate;  therefore, 
Kv  =  0. 

•  Water  in  the  reservoir  is  maintained  at  a  constant  level  using  stack 
coolant  water  recirculation. 

•  Water  reservoir  and  pump  outlet  are  opened  to  atmosphere; 
therefore,  inlet  pressure  equals  outlet  pressure  (Pin  =  Pout). 

•  Pressure  drop  across  the  stack  is  negligible. 

The  resulting  form  of  Eq.  (5)  is 

P,=pg(z2-z,)-KpQ2  (7) 

It  is  noted  that  Eq.  (7)  for  the  load  pressure  (Pi)  contains  two 
terms,  namely:  (i)  pressure  head  (g(z2  -  Zi))  due  to  height  (I);  (ii) 
the  pressure  head  due  to  pipe  roughness  ( I<PQ 2)  of  length  L. 

Now,  dividing  Eqs.  (1)  and  (2)  by JCom  and  Jw,  respectively,  and 
then  substituting  Eqs.  (3),  (4),  (6)  and  (7)  for  Tp,  Hp,  re  and  P\  in 
Eqs.(l)  and  (2),  yields 

CltWr  =  CW - Tf  +B(Mr  +  61Q2  +b2Q_(0r  +  b3C02  (8) 

Of  jcom 

=  a 1 Q2  +  a2Qcor  +  a3co 2  -  -  Zl  (9) 


Table  1 

Parameters  of  Tuscan  motor  and  Grundfos  CR2-30  pump 


a  Calculated  for  the  casing  volume  and  radius  of  150  cm3  and  4  cm,  respectively. 


Table  2 

Comparison  between  experimental  and  predicted  values  of  Grundfos  CR2-30  motor- 
pump  characteristic  data 


Current  (A) 

Experimental  values 

Predicted  values 

Speed  (rps) 

Flow  rate  (lps) 

Speed  (rps) 

Flow  rate  (lps) 

2.58 

214 

0.375 

211 

0.478 

3.47 

230 

0.497 

229 

0.545 

5.25 

262 

0.749 

261 

0.715 

6.42 

279 

0.854 

279 

0.785 

7.55 

295 

0.973 

295 

0.865 

where  B  —  —  B/Jcom.  hi  —  hi //com >  b2  —  — b2//Com»  b3  —  —b3/JCom, 
a]  =  Kp  -  ai/Jw,  a2  =  a2//w  and  a3  =  a3//w  are  the  modified  coef¬ 
ficients;  (z2  -  Zi)  is  the  distance  (height)  between  the  stack  inlet 
and  water  surface  of  the  reservoir  and  is  assumed  as  1  m. 

The  constants  C,  rf  and  Jcom  of  this  motor-pump  are  given  in 
Table  1.  The  moment  of  inertia1  of  water  yw)  is  calculated  for  a 
casing  volume  of  150  cm3  (liquid  volume  =  casing  volume)  and  a 
radius  of  impeller  of  4  cm,  respectively.  The  modified  coefficients 
in  Eqs.  (8)  and  (9)  are  obtained  by  simultaneously  solving  them  at 
steady  state.  For  this  purpose,  the  input-output  data  of  a  Tuscan 
36  V  motor  Grundfos  CR2-30  centrifugal  pump  given  in  Table  2  are 
utilized.  Five  input-output  data  points  are  taken  at  different  motor- 
pump  speeds  by  discretizing  the  characteristic  curves  given  in  [20]. 
The  estimated  values  of  coefficients  in  Eqs.  (8)  and  (9)  are 

[B  5t  b2  b3]  =  [0.09638  3.56801  0.02372  - 

0.0006531]  x  103 

[dj  d2  a3]  =  [8.295961  -0.060204  0.0001167]  x  105 

Simulation  of  the  dynamic  motor-pump  model  is  performed 
using  MATLAB-Simulink  with  the  motor-armature  current  as  input 
to  the  model.  The  motor-shaft  speed  and  pump  flow  rate  values 
obtained  from  the  model  are  listed  in  Table  2  for  comparison  with 
the  experimental  values  [20].  The  results  are  in  good  agreement 
with  the  experimental  values  of  shaft  speed  and  pump  flow  rate. 
Later,  this  dynamic  model  along  with  the  proposed  stack  thermal 
model  is  integrated  with  the  PEMFC  stack  model  for  control  stud¬ 
ies.  The  thermal  management  model  is  described  in  the  following 
section. 


3.  Thermal  management  model 

A  lumped  parameter  control-orientated  thermal  management 
model  of  a  PEMFC  stack  is  developed  based  on  the  energy  balance 
discussed  in  Section  1.  In  this  model,  the  heat  loss  to  the  surround¬ 
ings  is  taken  into  account  explicitly  rather  than  being  evaluated  in 
an  implicit  form  [16].  The  developed  thermal  model  can  predict 
the  temperature  of  the  stack,  the  reactants  at  exit  and  the  outlet 
coolant  water  for  a  given  set  of  load  (current),  inlet  coolant  water 
and  reactant  flows. 


1  Mass  x  radius2. 
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3.1  Control-orientated  stack  thermal  model 

The  development  of  the  thermal  model  is  based  on  physico¬ 
chemical  knowledge  of  the  fuel  cell  stack  and  with  following 

assumptions: 

(1)  the  temperature  is  uniform  over  the  whole  length  of  the  stack 
due  to  high  thermal  conductivity  and  the  provision  of  a  suffi¬ 
cient  number  of  cooling  plates  in  the  stack; 

(2)  pure  H2  gas  enters  at  the  anode  side  of  the  stack; 

(3)  no  liquid  water  enters  at  the  anode  and  cathode  sides; 

(4)  the  formation  of  product  water  at  the  cathode  is  in  vapour 
phase; 

(5)  water  transport  across  the  membrane  is  in  the  vapour  phase 
only. 


In  the  case  of  the  dead  end  on  the  anode  side  as  stated  in  Section 
3.2,  the  species  existing  at  anode  side  are  zero  despite  the  regular 
purging  of  hydrogen  at  the  anode  exit. 

The  sensible  heat  at  the  cathode  side  is  considered  for  four 
species,  namely:  oxygen,  nitrogen,  water  vapour  and  liquid  water 
at  the  cathode  exit.  From  assumptions  (4)  and  (5),  since  the  prod¬ 
uct  water  and  the  net  water  transport  across  the  membrane  are 
assumed  to  be  in  the  vapour  phase,  the  partial  condensation  of 
water  vapour  takes  place  on  the  cathode  side.  The  condensation 
occurs  when  the  partial  pressure  of  water  vapour  reaches  the  sat¬ 
uration  pressure  at  the  stack  temperature.  The  maximum  water 
vapour  that  is  held  by  the  gas  mixture  is  calculated  as  follows: 

P-ut 

Nw,v,c,out,max  —  (No2,c  ,out  +  N  n2  ,c,out )  “  pSat  > 

1  c’out  ~~  1  rc,out 


From  the  law  of  conservation  of  energy,  the  chemical  energy  of 
reactants  FI2-  02  equals  the  sum  of  the  electrical  energy  of  the  fuel 
cell  stack  and  the  thermal  energy  generated  due  to  overvoltage.  This 
thermal  energy  is  dispensed  in  the  form  of  latent  heat  of  vapouriza¬ 
tion  of  water  at  the  cathode  side,  heat  loss  to  the  surroundings  from 
the  stack,  and  sensible  heat  of  reactant  gases  and  coolant  water.  The 
energy  balance  is  given  by 

Qtheo  —  Qelec  +  Qsens  +  Qlatent  +  Qloss  (10) 

The  theoretical  energy  from  the  electrochemical  reaction  in  a 
PEMFC  is  calculated  by  the  product  of  the  heat  of  reaction  and  the 
number  of  moles  of  the  hydrogen  consumed  [21].  The  associated 
equation  is 

Q:heo  =  ^H2iCons  ^  rxn  (11) 

The  electrical  energy  produced  from  the  PEMFC  stack  with  ‘n’ 
single  cells  is  calculated  as 

Qeiec  =  nWst  (12) 

The  energy  loss  to  the  surroundings  due  to  natural  convection 
is  calculated  as 


Nw,c  ,out  =  NW(C  ,in  +  Ntrans  +  Nw.prod 

IfNw.c  ,out  ^  l^W,v,c,out,max  •  ^W.l.out  —  Nw, out  —  Nw,v,c,out,max  * 
Nw,v,c,out  —  ^W,c,out,max 

IfNw.c, out  —  Nw,v,c,out,max  •  ^W.l.out  =  0,  Nw,v,c,out  —  ^W.c.out 

The  above  equations  along  with  the  assumption  that  no  liquid 
water  enters  at  the  cathode  are  used  to  write  the  sensible  heat  on 
cathode  side  as 

Qsens, c  =  No2,c,outCp,  o2(Tc  ,out  —  Tq)  +  Nw,v,c,outCp?H20,v  (Tc  ,out  —  To) 

+  ^W,l,c,out^p,H20,  \(TC  ,out  —  T0)  +  JVn2  ,C,OUt  Cp,  N2{Tc  ,out  —  To) 
~  No2,c,inCp,02  (Tc  in  —  Tq)  —  Nw,v,c,inCp,H20,v(TClin  —  Tq) 

—  ^N2,c,in^p,N2  (Ic.in  ~  To)  (17) 

Finally,  the  total  sensible  heat  is  computed  as 

Qsens  —  Qsens, a  +  Qsens, c  +  Qsens, W  (18) 


Qloss 


—  hstAst{Tst  —  Tatm) 


(13) 


where  an  explicit  form  is  used  in  contrast  to  the  implicit  form  that 
was  provided  in  [  16]  to  account  for  the  heat  loss  to  the  surroundings. 
The  heat  loss  in  the  above  equation  is  calculated  using  the  value  of 
the  film  heat  transfer  coefficient  (hst)  obtained  from  the  Nusselt 
number.  The  Nusselt  number  for  a  vertical  plate  has  been  derived 
by  Churchil  and  Chu  [22,23],  i.e. 


0.825  +  0.387- 


fia01 


[1  +  (0.492/Npr) 


,0.56,8/27 


(14) 


from  which  hst  is  obtained. 

The  sensible  heat  of  coolant  water  stream  is  calculated  using  the 
following  formula  where  the  loss  of  liquid  mass  by  vapourisation 
is  assumed  to  be  negligible  at  the  exit  of  the  coolant  water  stream: 


Qsens, W  =  NWCPiw(lw,out  -  Tw in) 


(15) 


From  assumptions  (2)  and  (3),  the  sensible  heat  at  anode  side 
is  considered  for  two  species,  namely  pure  hydrogen  and  water 
vapour.  The  amount  of  liquid  water  on  the  anode  side  is  neglected 
by  virtue  of  a  predominant  electro-osmotic  drag.  Under  these 
assumptions,  the  energy  balance  on  anode  is  written  as 


Qsens, a  —  l^H2,a,outCp,H2  (la, out  —  Tq)  +  Nw,v,a,outQj,H20,v(la,out  —  Tq) 


The  latent  heat  of  water  vapour  due  to  phase  change  is  calcu¬ 
lated  by  noting  that  the  molar  flow  rate  of  water  involved  during 
the  phase  change  is  the  difference  between  the  total  water  on  the 
cathode  side  and  the  maximum  water  vapour  that  is  held  by  the 
gas  mixture  (Nw,v,c,out  ,max)* 

If  Nw.c.out  >  Nw,v,c  , out, max  • 

Qlatent  —  (Nw,c,out  —  ^W,v,c, out, max  ^vapourisation  (19) 

IfNw.c, out  —  l^W,v,c,out,max  • 

Qlatent  =  9  (20) 

where  Hvapourisation  is  the  latent  heat  of  vapourization  of  water  given 
by 


^vapourisation  =  45070  -  41 .9 Tst  +  3.44  x  lO"3^  +  2.54  X  10-6T3 
-  8.98  X  10-10rs4t 


3.11  Calculation  of  stack  temperature 

The  overall  energy  balance  equation  of  fuel  cell  stack  is  given  by 
Energy  accumulation  =  Energy  in  -  Energy  out 

MstCptst-^jr  —  Qtheo  —  (Qelec  +  Qsens  +  latent  +  Qloss)  (21 ) 


^H2,a,inQj,H2(la,in  Tq)  ^W,v,a,in^p,H20,v(la,in  Tq) 


(16) 


Eq.  (21 )  is  vital  to  the  calculation  of  the  dynamic  stack  tempera¬ 
ture  where  Eqs.  (11 ),  (12),  (13),  (18)  and  (19)  are  invoked  to  compute 
Qtheo*  Qelec*  Qloss*  Qsens  and  Qlatent.  respectively. 
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3.1.2.  Calculation  of  coolant  water  outlet  temperature 

The  outlet  temperature  of  the  coolant  water  (Tw.out)  is  calculated 
by  equating  the  sensible  heat  of  coolant  water  to  the  convective  heat 
transfer  of  coolant  water,  i.e. 


^wCP)h2o( ^w.out  -  Tw.in)  =  ^w^wC^lmtd)  (22) 

where  the  logarithmic  mean  temperature  difference  (Tlmtd) 
between  cooling  water  and  stack  is  used  to  represent  the  average 
temperature  across  the  inlet  and  outlet  of  the  stack  and  is  given 
by 


Tlmtd  = 


Tw.out  T^v.in 


ln((Tst  —  Tvvtin)/(Tst  —  Tw.out)) 


(23) 


After  substituting  for  TLM td  and  by  rearranging  Eq.  (22),  the  final 
expression  for  coolant  outlet  temperature  takes  the  form: 


Tw.out  —  Tst  — 


ln(Tst  -  Tvv.in) 


hw^w 

NwCp,H2o 


(24) 


3.1.3.  Calculation  of  cathode  and  anode  exit  gas  temperatures 

The  temperatures  of  the  exit  gas  flow  streams  at  the  anode  and 
the  cathode  are  calculated  by  equating  the  sensible  heat  of  each 
gas  mixture  to  the  convective  heat  transfer  of  the  gas  mixture  on 
each  side  (anode,  cathode).  The  average  exit  temperatures  of  the  gas 
mixture  at  the  anode  and  the  cathode,  Ta(0ut  and  TCj0Ut,  are  obtained 
using  the  expressions  reported  in  [16]: 


Ta.out  —  2 


Tst- 


Qsens.a  +  Qma 


(hA)  a 


(25) 


Tc.out  —  2 


Tst- 


Qsens.c  +  Qlatent.c  —  Qma: 


(hA)c 


(26) 


where  Qmass.a  and  Qmass.c  are  the  energy  changes  due  to  mass  trans¬ 
fer  and  consumption  of  reactants  at  the  anode  and  the  cathode  side, 
respectively,  namely 


Qmass.a  —  ^transCp?H20.v(T'st  —  Tq)  +  NH2,ConCp,H2(T'st  —  Tq), 
Qmass.c  =  NtransCp?H20.v(T"st  —  To)  +  NH2,conCp,H20,l(T"st  —  To) 
—No2<conCp,02(.Tst  —  To) 


This  completes  the  development  of  PEMFC  stack  thermal  man¬ 
agement  model.  The  temperatures  of  the  stack,  the  coolant  outlet 
water  and  the  anode  and  cathode  exit  streams  can  be  obtained  from 


Eqs.  (21),  (24),  (25)  and  (26),  respectively.  The  developed  thermal 
model  is  now  integrated  with  the  PEMFC  stack  and  air  flow  com¬ 
pressor  models  coded  in  MATLAB-Simulink  environment  [24]  to 
carry  out  a  study  of  the  thermal  effects  on  the  PEMFC  system.  Due 
to  the  result  of  thermal  model  integration,  an  input  (motor-pump 
current)  and  output  (stack  temperature)  are  added  in  addition  to 
the  input  (compressor  motor  voltage)  and  output  (air  flow  rate) 
of  the  PEMFC  stack  and  air  flow  compressor  models.  The  equa¬ 
tions  describing  the  PEMFC  stack  and  the  air  flow  compressor  are 
provided  in  Appendix  A.  A  schematic  of  the  integrated  model  is 
shown  in  Fig.  1 .  Further,  this  integrated  model  is  used  for  interac¬ 
tion  analysis,  as  discussed  in  Section  4.4.  An  overview  of  the  PEMFC 
system  with  which  the  thermal  model  is  integrated  is  given  below 
for  convenience. 

3.2.  Review  of  PEMFC  stack  model 

Pukrushpan  et  al.  [17]  have  developed  a  system  level  lumped 
parameter  model  based  on  physical-chemical  knowledge  of  the 
processes  involved  in  the  PEMFC  stack  and  its  reactant  (air)  sup¬ 
ply  system  [17,25].  It  contains  the  air  flow  compressor  model  and 
the  semi-empirical  current-voltage  model  (i.e.  polarization  curve). 
This  model  is  employed  to  maintain  the  oxygen  excess  ratio  at  the 
desired  value  of  2  when  the  stack  is  subjected  to  dynamic  load 
changes.  The  salient  features  and  the  underlying  assumptions  of 
this  model  are  as  follows: 

•  The  parameters  used  in  PEMFC  stack  model  are  based  on  the  75 
kW  stack  used  in  the  FORD  P2000  fuel  cell  prototype  vehicle. 

•  A  semi-empirical  model  of  the  current-voltage  curve  is  built 
based  on  a  wide  range  of  stack  operating  temperature  (50-90  °C) 
data. 

•  The  stack  comprises  381  cells  of  active  area  280  cm2. 

•  The  humidity  of  the  reactant  gases  is  constant. 

•  The  fuel  (H2)  supply  on  anode  side  is  considered  to  be  static  and 
only  the  air  flow  rate  dynamics  are  considered. 

•  Hydrogen  is  supplied  at  the  anode  side  in  stoichiometry  propor¬ 
tion  and  is  adjusted  instantaneously  with  a  proportional  control 
valve  placed  at  the  high-pressure  cylinder. 

•  The  dead  end  is  considered  at  the  anode  exit. 

•  The  temperature  of  the  stack  is  assumed  to  be  constant  at  80 0  C 
and  is  uniform  over  the  whole  length  of  the  stack. 

•  All  gases  are  assumed  to  possess  ideal  gas  behaviour. 

•  The  temperature  on  the  cathode  flow  channel  is  assumed  to  be 
equal  to  the  stack  temperature. 
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Table  3 

Parameters  to  be  supplied  to  simulate  integrated  stack  model 


a  Area  of  anode,  cathode  and  single  cooling  plate  are  calculated  based  on  Marl<902 
Ballard  stack  data. 

3.3.  Simulation  of  integrated  model 

Additional  parameters  introduced  due  to  the  integration  of  the 
thermal  model  are  set  to  the  values  given  in  Table  3.  The  settings  for 
the  air  flow  compressor  and  voltage-current  model  are  set  to  their 
suggested  values  [17].  The  inputs  required  to  simulate  the  stack 
thermal  model  are  the  molar  flow  rates  of  hydrogen,  oxygen,  nitro¬ 
gen,  water  streams,  water  transport  across  the  membrane  and  total 
pressure  at  the  cathode  exit,  and  electrical  energy  generated  from 
stack.  These  inputs  are  taken  from  the  anode,  the  cathode,  the  mem¬ 
brane  hydration  and  current-voltage  models  of  the  PEMFC  stack 
and  air  flow  compressor  model. 

4.  Results  and  discussion 

The  developed  thermal  model  is  validated  based  on  the  results  of 
Ballard  V  5  kW  thermal  data  available  in  the  literature  [26].  Due  to 
the  non-disclosure  of  Ballard  900  stack2  data,  the  input  conditions 
of  5  kW  Ballard  stack  are  scaled-up  to  facilitate  the  comparison  with 
the  results  from  the  simulation  of  the  70  kW  stack  under  consider¬ 
ation.  The  scaling  is  done  by  maintaining  similar  temperatures  of 
coolant  water  and  inlet  feeds  at  the  anode  and  the  cathode.  Table  4 
provides  this  information  along  with  the  values  of  the  physical  data 
used  for  validating  the  model.  Further,  the  stack  heat  capacity,  oxy¬ 
gen  and  inlet  coolant  water  flow  rates  of  the  5  kW  unit  [26]  are 
scaled  by  factor  of  14.  The  inlet  hydrogen  flow  rate  at  the  anode  is 
taken  as  per  stoichiometry.  The  aforementioned  variables  are  lin¬ 
early  scaled  based  on  the  assumption  that  the  polarization  curves 
of  both  Ballard  V  (5  kW)  and  Ballard  900  (70  kW)  stacks  are  similar. 
This  assumption  is  valid  since  the  same  membrane  (Nafion-117)  is 
used  in  both  the  stacks. 

The  stack  and  exit  temperatures  of  the  outlet  streams  of  the 
70  kW  stack  would  be  equivalent  to  those  of  the  5  kW  stack,  pro¬ 
vided  the  thermal  energy  released  from  70  kW  stack  equals  14 


2  Peak  power  of  the  stack  is  70  kW  and  the  latter  contains  381  cells  of  active  area 
280  cm2. 


Table  4 

Table  of  Ballard  V  (5  kW)  and  scaled  data  of  70  kW  (Ballard  900)  for  similar  inputs 


times  the  thermal  energy  released  from  the  5  kW  stack.  The  amount 
of  load  required  to  generate  this  thermal  energy  is  calculated 
as 

^chemical  —  ^electrical  =  14  x  (Fchemical  —  ^electrical)  (27) 

381  x  1st  (^  -  0.82)  =  14  x  35  x  20  (|^  -  O.82)  (28) 

Further,  the  parameters  ((M)a,  (M)c,  (ftA)w)  of  70  kW  stack  are 
obtained  by  linearly  scaling  the  Ballard  V  (5  kW)  stack  values.  The 
heat  capacity  (MstCp>st)  of  the  70  kW  stack  is  taken  as  14  times  that 
of  the  5  kW  stack  to  facilitate  the  same  settling  time  as  that  of  the 
5  kW  stack  temperature.  To  validate  the  developed  model  with  the 
experimental  values,  the  data  given  in  Table  4  is  used  instead  of  the 
corresponding  same  parameters  given  in  Table  3. 

The  heat  transfer  coefficient  (hst)  between  the  air  and  the  stack 
surfaces  is  calculated  using  Eq.  (14)  and  is  shown  in  Fig.  2.  The 
resulting  values  are  found  to  be  in  the  range  of  the  experimental 
values  (3-7  W  m_2K_1)  given  in  [27]  when  the  stack  is  oper¬ 
ated  between  50  and  90 0  C.  The  steady-state  temperatures  of  the 
stack,  coolant  water  and  cathode  outlet  of  the  70  kW  stack  obtained 
from  the  simulation  of  the  developed  thermal  model  are  shown  in 
Table  5  for  the  same  inlet  temperatures  of  the  coolant  water  and 
the  reactants  as  that  of  the  5  kW  stack.  The  last  column  of  this 
table  contains  the  corresponding  outputs  of  the  5  kW  for  the  pur- 


Fig.  2.  Free  convective  heat  transfer  coefficient  from  stack  to  surroundings. 
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Table  5 

Validation  of  thermal  model  with  Ballard  V  (5  kW)  temperature  data 


a  Settling  time  of  stack  temperature. 


pose  of  comparison.  It  is  seen  that  the  predicted  values  match  well 
with  the  experimental  data.  The  transient  temperature  responses  of 
the  stack,  cathode  and  coolant  outlet  streams  are  shown  in  Fig.  3. 
It  is  inferred  from  Fig.  3  that  the  settling  time  of  the  stack  tem¬ 
perature  is  60  min  while  that  of  the  Ballard  V  stack  is  55  min. 
The  discrepancy  between  experimental  and  predicted  tempera¬ 
tures  and  the  settling  time  is  due  to  the  difference  in  cell  voltages 
of  the  5  and  70  kW  stacks  which  is  assumed  to  be  constant  during 
the  load  calculation.  The  temperature  of  the  anode  exit  gas  is  not 
reported  here  due  to  the  dead-end  configuration  of  the  70  kW  stack, 
as  stated  in  Section  3.2,  which  is  open-end  in  the  5  kW  Ballard  V 
stack. 

The  results  of  the  pump  model  validation  were  presented  in 
Section  2.  Thus  the  proposed  model  provides  predictions  that  are 
concurrent  with  the  physics  of  the  process  as  reported  in  the  lit¬ 
erature.  The  thermal  model,  water  pump  model  and  the  existing 
PEMFC  stack  model  are  integrated,  as  described  in  Section  3.3.  This 
integrated  model  is  used  to  examine  the  effect  of  (i)  the  number  of 
cooling  plates  on  the  coolant  outlet  temperature  and  (ii)  the  stack 
heat  capacity  on  the  settling  time  of  the  stack  temperature.  The 
load  current  and  compressor  motor  voltage  are  set  to  191 A  and 
164  V,  respectively,  and  result  in  a  net  power3  of  40 kW  [25].  The 
stack  is  operated  at  353.15  K.  The  motor-armature  current  of  the 
water  pump  is  set  to  1.157  A  to  ensure  that  the  stack  temperature  is 
maintained  at  353.15  K.  The  parameters  given  in  Table  3  are  used  in 
conjunction  with  these  operating  conditions.  Further,  the  effect  of 
stack  temperature  on  the  settling  times  of  air  flow,  the  stack  volt¬ 
age,  the  net  power  and  the  oxygen  excess  ratio  under  closed-loop 
operation  is  studied. 


Number  of  cooling  plates 
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Fig.  4.  Effect  of  number  of  cooling  plates  on  exit  coolant  water  temperature. 


4.1.  Influence  of  number  of  cooling  plates  on  exit  coolant  water 
temperature 

The  main  advantage  of  this  study  is  that  the  parasitic  losses  due 
to  auxiliary  equipment  like  radiators  can  be  reduced.  These  radia¬ 
tors  are  used  to  cool  the  outlet  coolant  water  to  a  reasonable  low 
temperature  before  it  is  fed  back  at  the  stack  inlet. 

The  number  of  cooling  plates  that  are  chosen  for  this  study  are 
in  the  range  of  n  /4  to  3 n  / 4,  where  ‘n’  is  the  number  of  single  cells 
in  the  stack.  From  Fig.  4  it  is  observed  that  the  temperature  of  the 
coolant  outlet  reaches  close  to  that  of  the  stack,  i.e.  the  temperature 
difference  between  the  stack  and  the  coolant  outlet  water  is  1 0  C 
or  below  as  the  number  of  cooling  plates  is  increased  beyond  220. 
A  further  increase  in  the  number  of  cooling  plates  leads  to  only  a 
marginal  decrease  in  the  temperature  difference  at  the  expense  of 
increase  in  weight  and  size  of  the  system.  In  view  of  these  observa¬ 
tions,  the  number  of  cooling  plates  is  fixed  at  220  for  subsequent 
studies. 

4.2.  Influence  of  stack  heat  capacity  on  settling  time  of  stack 
temperature 

Eq.  (21 )  describes  the  dynamic  relationship  of  stack  temperature 
with  Mst  and  Cp,st.  Simulation  studies  are  conducted  to  quantify  the 
influence  of  mass  (Mst)  and  average  specific  heat  capacity  (Cp,st)  of 
the  stack  on  the  settling  time  of  the  stack  temperature  which,  in 
turn,  effects  the  stack  start-up  time. 

The  settling  time  of  the  stack  temperature  is  studied  by  grad¬ 
ually  decreasing  the  heat  capacity  quantity,  i.e.  (MCp)st,  below 
the  nominal  value  of  the  70  kW  stack.4  The  heat  capacity  can 
be  varied  by  choosing  materials  alternative  to  graphite  for  bipo¬ 
lar  or  cooling  plates  that  result  in  lower  values  of  (MCp)st.  From 
Fig.  5,  it  is  observed  that  the  settling  time  of  the  stack  temperature 
decreases  almost  linearly  with  decrease  in  the  stack  heat  capacity, 
as  expected.  Based  on  this  observation,  fuel  cells  of  today  are  built 
using  a  lower  heat  capacity  than  graphite.  Coated  bipolar  plates 
and  cooling  flow-field  plates  made  out  of  aluminum  possess  higher 
thermal  and  electrical  conductivity  than  graphite  in  addition  to  a 
low  heat  capacity.  Therefore,  aluminum  can  be  preferred  for  making 
bipolar  (coated)  and  cooling  plates.  This  study  has  been  carried  out 


3  Net  power  =  stack  power  -  power  drawn  by  the  air  compressor. 


4  The  product  of  (MCp)st  is  taken  as  490  kj  K-1  based  on  the  Ballard  5  kW  stack. 
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Fig.  5.  Effect  of  heat  capacity  of  stack  on  settling  time  of  stack  temperature. 


on  the  lines  of  quantifying  the  influence  of  the  material  on  the  set¬ 
tling  times  of  stack  temperature  as  this  which  is  important  during 
start-up  of  the  stack. 

4.3.  Effect  of  stack  temperature  on  air  flow  compressor,  stack 
voltage,  oxygen  excess  ratio  and  net  power  transients 

To  investigate  these  effects,  studies  are  conducted  under  two 
different  conditions,  namely:  (i)  constant  temperature  (as  per  the 
assumption  made  in  [25]  and  re-stated  in  Section  3.2)  and  (ii)  tem¬ 
perature  is  controlled  at  a  desired  value  (353.15  I<).  The  difference 
between  these  two  conditions  is  that  in  the  second  case,  the  stack 
temperature  and  the  controller  dynamics  are  taken  into  account. 
The  results  obtained  for  both  the  cases  are  discussed  below. 

The  abovementioned  cases  are  studied  in  a  closed  loop  with  the 
same  load  (current)  profile  as  given  in  [17]  and  shown  in  Fig.  6.  In 
the  first  case  (constant  temperature),  a  PI  controller  is  employed  to 
control  the  air  flow  rate  such  that  an  oxygen  excess  ratio  of  two  is 
maintained.  The  integrated  model  is  used  to  simulate  the  second 
case  where  temperature  dynamics  are  taken  into  consideration. 
Here,  two  single-loop  PI  controllers  are  used  to  control  the  air  flow 
rate  and  stack  temperature  separately  by  manipulating  the  com¬ 


Time  (sec) 

Fig.  7.  Dynamic  response  of  air  flow  rate  under  various  load  conditions. 

pressor  motor  voltage  and  motor-armature  current,  respectively. 
The  resulting  settling  times  of  outputs  (voltage,  air  flow  rate)  and 
performance  variables  (oxygen  excess  ratio,  net  power)  are  com¬ 
pared  for  both  cases. 

The  dynamic  response  of  the  air  flow  rate,  oxygen  excess  ratio, 
stack  voltage,  net  power  and  stack  temperature  are  shown  in 
Figs.  7-11  ,  respectively.  It  is  observed  from  Figs.  7  and  8  that  the 
air  flow  rate  and  oxygen  excess  ratio  settle  within  4  s  in  both  cases. 
On  the  other  hand,  it  can  be  noticed  from  Fig.  9  that  the  stack  volt¬ 
age  shows  a  significant  bias  when  compared  with  the  results  of  the 
constant-temperature  case  where  it  settles  within  5  s.  The  shift 
in  voltage  away  from  the  steady-state  value  due  to  temperature 
change  can  be  expected  to  result  in  a  bias  in  the  stack  net  power 
when  compared  with  the  constant  temperature  case,  as  shown  in 
Fig.  10.  The  difference  in  the  net  power  output  is  in  the  range  of 
250-300  W,  which  may  not  be  discernible  due  to  the  scaling  that  is 
used  to  display  the  net  power  (in  kW).  The  bias  in  voltage  is  due  to 
the  fact  that  the  voltage  is  a  function  of  stack  temperature  which 
does  not  settle  within  the  duration  spanned  by  two  successive  load 
changes,  the  stack  temperature  and  the  corresponding  input  motor- 
armature  current  profiles  are  shown  in  Fig.  11.  This  is  associated 
with  the  fact  that  the  open-loop  settling  time  of  the  stack  temper- 
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Fig.  12.  Closed-loop  settling  time  of  stack  temperature  for  step-change  in  load  from 
191  to  230  K. 


Fig.  9.  Dynamic  response  of  stack  voltage  under  various  load  conditions. 


ature  is  about  50  min  (3000  s)  and  therefore  it  is  only  reasonable 
to  design  a  controller  that  can  result  in  closed-loop  settling  times 
greater  than  5  min  (1  /10th  of  open-loop  settling  time).  In  fact,  with 
the  tunings  used  in  the  simulations  the  settling  time  of  the  stack 
temperature  is  6  min  (nearly  300  s)  as  shown  in  Fig.  12  for  a  step 
change  in  load  (191-230  A). 

Since  the  air  flow  is  one  of  the  carriers  of  thermal  energy  from 
the  stack,  it  influences  the  stack  temperature.  By  contrast,  the  effect 
of  the  stack  temperature  on  the  compressor  air  flow  is  virtually 
non-existent.  In  order  to  quantify  the  extent  of  interactions  among 
the  input-output  pairs,  interaction  studies  are  conducted  using  the 
Relative  Gain  Array  (RGA)  tool.  The  results  are  presented  in  the 
following  section. 

4.4.  Interaction  studies— from  a  control  perspective 

The  primary  objective  of  this  study  is  to  assess  the  extent 
of  interaction  for  control  purposes.  The  steady-state  interaction 
between  inputs  (compressor  motor  voltage  and  motor-pump  cur¬ 
rent)  and  outputs  (air  flow  rate  and  stack  temperature)  of  interest 
are  studied  using  the  integrated  model  built  in  a  MATLAB-Simulink 
environment.  The  interaction  between  each  input-output  pair  is 
facilitated  by  the  well-known  RGA  analysis  [28],  i.e. 


Fig.  11.  Profiles  of  motor  armature  current  (manipulated)  and  stack  temperature 
(controlled)  when  controlled  at  353.15  K. 


RGA  = 


ux  u2 

1.078  -0.074 

-0.079  1.062 


Vi 

V2 


The  RGA  matrix  suggests  that  the  pairs  output  1  (y i :  air  flow 
rate)-input  1  (ux:  compressor  motor  voltage)  and  output  2  (y2: 
stack  temperature )-input  2  (u2:  motor-pump  current)  have  a  strong 
interaction.  On  the  other  hand,  the  off-diagonal  pairs,  namely,  input 
1 -output  2  and  input  2-output  1  exhibit  very  weak  interaction. 
Therefore,  it  is  concluded  that  the  loops  can  be  controlled  on  an 
individual  basis. 

5.  Conclusions 

In  this  investigation,  a  control-orientated  thermal  management 
model  together  with  a  dynamic  water  pump  model  has  been  devel¬ 
oped  to  predict  the  temperatures  of  the  stack,  cathode,  anode  and 
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coolant  exit  streams  for  changes  in  air  flow  compressor  voltage, 
coolant  water  flow  rate  (motor-pump  armature  current)  and  load 
conditions.  The  practical  utility  of  these  models  has  been  demon¬ 
strated  by  integrating  them  with  PEMFC  stack  current-voltage  and 
air  flow  compressor  models  that  exist  in  the  literature. 

The  results  obtained  from  the  simulation  studies  of  the  inte¬ 
grated  system  showed  good  agreement  with  the  experimental  data 
available  in  the  literature.  Further,  these  results  were  shown  to  be 
useful  in  obtaining  good  physical  insights  into  (i)  the  choice  of  alter¬ 
native  materials  for  cooling  plates  and  bipolar  plates  other  than 
graphite  that  can  reduce  the  settling  time  of  the  stack  temperature 
and  the  overall  size  of  the  PEMFC  system,  and  (ii)  the  required  num¬ 
ber  of  cooling  plates  to  ensure  that  the  temperature  of  the  coolant 
outlet  water  is  close  to  that  of  the  stack. 

The  open-loop  study  of  the  thermal  subsystem  showed  that  the 
settling  time  of  the  temperature  is  about  50  min,  for  a  step-change 
in  the  load.  The  closed-loop  system  is  tuned  such  that  the  stack 
temperature  settles  in  about  6  min.  Interaction  studies  were  car¬ 
ried  out  to  analyze  the  effect  of  the  thermal  sub-system  on  the  air 
flow  rate  sub-system  and  vice  versa.  These  studies  revealed  that  the 
settling  times  of  most  of  the  variables  remain  unaffected  due  to  the 
dynamics  of  the  temperature  loop,  while  the  voltage  and  net  power 
settle  only  after  the  temperature  settles  down.  The  implication  of 
this  fact  from  a  control  perspective  is  that  the  two  sub-systems 
can  be  treated  as  non-interacting  systems  suggestive  of  the  use  of 
single-loop  controllers  for  each  of  these  loops. 

The  integrated  model  developed  in  this  paper  is  useful  in  a 
few  important  aspects.  For  instance,  it  can  be  employed  to  study 
process  heat  integration  and  disturbance  (load)  rejection  to  avoid 
oxygen  starvation  and  system  faults  detection  and  diagnosis.  Fur¬ 
ther,  various  control  algorithms  can  be  tested  using  this  model  for 
maintaining  the  stack  temperature.  Advanced  control  schemes  may 
be  required  to  minimize  the  settling  times  of  the  air  flow  rate  and 
to  minimize  the  fluctuations  in  the  stack  temperature.  These  issues 
provide  directions  for  future  studies. 


Appendix  A 

A  A.  Brief  overview  of  compressor  and  polarization  curve  model 


Supply  manifold.  The  inlet  mass  flow  is  compressor  air  flow  Wcp 
and  the  outlet  mass  flow  is  Wsmiout.  State  equations  are 

=  Wcp  -  Wsm,out  (A.4) 

. .  =  77 — ( WcpT'cp.out  —  VVsnpoutTsm)  (A. 5 ) 

at  vsm 

where  VSm,  Y,  Ra  and  Tsm  are  the  supply  manifold  volume,  specific 
heat,  ideal  gas  constant  and  supply  manifold  temperature  of  air, 
respectively. 

Humidifier  (static)  model.  Before  the  air  is  fed  into  the  humidifier, 
it  is  cooled  to  a  desired  temperature  in  the  static  cooler.  Wvcl  is  the 
amount  of  water  vapour  present  in  the  air-water  vapour  mixture 
before  the  air  is  humidified.  Wv  hm,  pvhm  and  </>h m  are  the  vapour 
flow  rate,  vapour  pressure  and  relative  humidity  of  the  gas  mixture 
at  the  humidifier  exit  as  follows: 


Wv,hm  —  Wv.cl  +  Wy  jnj 

(A.6) 

„  Wv hm  Ma  M 

Pv,hm  -  Wad  MvP a,cl 

(A.7) 

i  Pv,hm 

9hm  ~  Psat(Thm) 

(A.8) 

where  Wv  inj  is  the  amount  of  water  vapour  injected  in  the  humidi¬ 
fier,  Wa  d  the  dry  flow  rate  of  air,  pa  cl  the  partial  pressure  of  dry  air, 
Thm  the  temperature  of  humidifier,  and  Ma  and  Mv  are  the  molecular 
mass  of  air  and  water  vapour,  respectively. 

Similar  equations  are  used  on  hydrogen  side. 

Cathode  flow  model.  The  model  is  developed  using  the  mass 
conservation  principle  and  thermodynamic  and  psychometric 
properties  of  air.  The  flow  model  is  lumped  as  a  CSTR,  i.e.  the  vari¬ 
ables  at  cathode  exit,  namely,  temperature  Tca, out,  pressure  pca,out 
and  oxygen  mole  fraction  y02,ca, out  are  assumed  to  be  the  same  as 
the  they  are  in  the  cathode  flow  channel,  Tca,  pca  and  yo2,ca-  The 
mass  continuity  equations  are 

=  W02,ca,in  -  W02,ca,out  -  W02,reacted  (A.9) 


dmN2,ca 

dt 


—  ^N2,ca,in  —  ^N2,ca,out 


(A.10) 


The  PEMFC  stack  system  contains  the  models  of  the  fuel  cell 
stack  and  its  auxiliary  components  such  as  compressor  model,  man¬ 
ifold  model  and  humidification  model  [17,25]. 

Compressor  model.  The  inputs  to  the  model  include  inlet  air 
pressure  pcpin,  its  temperature  rcPiin,  the  voltage  command  to  the 
compressor  motor  vcm,  and  downstream  pressure  psm.  The  only 
dynamic  state  in  this  model  is  compressor  speed  oo cp.  The  torque 
required  to  drive  the  compressor  is  given  by 


Cp  Tatm 
&>cp  Pep 


Psm\{y-')/r 
Patm  / 


wcl 


(A.l ) 


The  dynamic  behaviour  of  the  compressor  speed  is  represented 
by 


Jcp- 


d  co. 


’cp 


dt 


l  —  Tcp) 


(A.2) 


dmw>ca 

dt 


—  Wv,ca,in  -  Wv.ca.out  +  Wv.ca.gen  +  tVv.membr  —  ^l.ca.out 


(A.ll) 


Return  manifold.  The  return  manifold  pressure  at  the  cathode  is 
modelled  by 

dprm  RaTnn,!.,  *A/  \ 

-t ~  =  —rj (Wcaiout  -  Wrm)0ut)  (A.l 2) 

dt  Vrm 

where  Vrm  is  the  return  manifold  volume  and  Trm  is  the  temperature 
of  the  gas  in  the  return  manifold. 

Anode  flow  model.  Assumptions  are  similar  to  the  cathode  flow 
models.  The  hydrogen  is  compressed  and  stored  in  a  hydrogen  tank 
and  the  anode  inlet  flow  rate  is  assumed  to  be  supplied  instan¬ 
taneously  with  a  valve  of  a  fast  actuator  to  maintain  a  minimum 
pressure  difference  between  the  anode  and  the  cathode.  The  state 
equations  are 


The  compressor  motor  torque  is  calculated  as 

tcm  =  P cm  „ — (vcm  —  /<v<^cp) 

Kcm 


=  WH2,an,in  -  W^.an.out  -  Wh2 .reacted 

(A.3)  dmw,a 


dt 


—  Wv,an,in  Wv.an.out  Wv.membr  ^1, an, out 


(A.13) 

(A.14) 


where  C,  Rcm  and  kv  are  motor  constants  and  rjcm  is  motor  mechan-  Membrane  hydration  model.  This  calculates  the  water  content  in 

ical  efficiency.  the  membrane  and  rate  of  mass  flow  across  the  membrane.  Both 
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water  content  and  mass  flow  are  assumed  to  be  uniform  over  the 
surface  area  of  the  membrane.  The  membrane  water  content  and 
the  rate  of  mass  flow  across  the  membrane  are  functions  of  the 
stack  current  and  relative  humidity  of  the  reactants  in  the  anode 
and  the  cathode  flow  channels. 

Water  transport  across  the  membrane  is  considered  to  occur 
in  two  ways,  namely:  (i)  electro-osmotic  drag  where  the  water 
molecules  are  dragged  from  the  anode  to  the  cathode  side  by 
protons;  (ii)  back  diffusion,  in  this  case  the  water  molecules  are 
transferred  due  to  the  concentration  gradient  of  water  across  the 
membrane.  The  expressions  are  as  follows: 

at  ^  at  n  t’v.ca  —  Cv,a 

F»v, osmotic  -^dri  —  —  ~ 

r  tm 

The  net  transport  of  water  across  the  membrane  is  given  by 


^v.membr  —  ^v,  osmotic  +  ^vdiff  (A.15) 

where  NVOSmotic  (mols-1  cm-2)  is  the  net  water  flow  from  anode 
to  cathode  of  one  cell  caused  by  electro-osmotic  drag,  i  (A cm-2) 
the  stack  current  density,  Fthe  Faraday  number  and  nd  the  electro- 
osmotic  drag  coefficient,  Nvd iff  (mols-1  cm-2)  the  net  water  flow 
from  cathode  to  anode  of  one  cell  caused  by  back-diffusion,  cv 
(mol  cm-3)  the  water  concentration,  tm  the  thickness  of  the  mem¬ 
brane  (cm)  and  Dw  (cm2  s-1 )  is  the  diffusion  coefficient  of  water  in 
the  membrane. 

Stack  voltage  model.  The  voltage  is  calculated  using  a  combina¬ 
tion  of  physical  and  empirical  relationships,  and  is  given  by 


Vfc  =  E  -  Vact  -  Vohm  -  Vconc  =  E  -  [v0  +  Va(l  -  e  c'')] 


[^ohm] 


i 


c3  "I 


(A.16) 


where  E  is  the  open-circuit  voltage,  vact.  v0hm  and  vCOnc  are  the 
voltage  drops  due  to  activation  loss,  ohmic  loss  and  concentration 
loss,  respectively.  The  open-circuit  voltage  is  calculated  using  the 
energy  balance  between  the  reactant  and  products,  and  the  Faraday 
constant.  It  is  given  by 


E  =  1.229-8.5  x  10-4(Tst- 298.15) +  4.3085 


X  10-5Tst[ln(pH2 )  +  0.5  ln(po2 )]  (A.17) 

where  v0  (V)  is  the  voltage  drop  at  zero  current  density,  i  the  current 
density  (A cm-2),  Rohm  the  internal  electrical  resistance  (£2 cm2), 
and  va,  /max,  Ci,  c2  and  c3  are  constants.  The  parameters  in  the 
expression  (A.16)  are  determined  using  nonlinear  regression  of  fuel 


cell  polarization  data  from  an  automotive  propulsion-sized  fuel  cell 
stack.  The  results  obtained  from  the  model  show  good  accuracy 
when  it  is  tested  at  an  operating  temperature  between  40  and  100 0 
C  and  for  various  cathode  pressures. 
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